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1 Introduction 

The physics of ultracold bosons confined in a double-well potential has attracted 
a great deal of attention since the theoretical prediction of Josephson-like oscilla- 
tions of the atoms population and the existence of self-trapped states '-2. In addi- 
tion, the recent experimental realization of a bosonic Josephson junction (BJJ) by 
the Heidelberg group using ^^Rb atoms^ has triggered the possibility of practical 
applications and extensions to other physical scenario s 

The theoretical prediction of Smerzi et. al. ' was made by means of the mean- 
field Gross-Pitaevskii (GP) equatio n ''^i^'i'^ , which correctly captures the tunnel- 
ing dynamics of the population and its coupling to the phase difference between 
the two sides of the barrier. A further simplification, which turned out to be partic- 
ularly useful, is the consideration of only the lowest two modes of the GP equation. 
Most of the semi-classical predictions of this two-mode approach '^''^'^, dealing 
with the Rabi to Josephson transition, have been confirmed in a Josephson exper- 
iment where the two modes are two distinct internal states of the atom'^. 
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The two-mode model can be requantized giving rise to a two-site Bose-Hubbard 
(BH) model^ii^^. It is worth noting that the regime of applicability of the quan- 
tized two mode approximation can extend further: recent examples are the exper- 
iments on BJJ, the production of number squeezed states, and a non-linear atom 
interferometer These phenomena are beyond GP, as they involve entangled 
states of the atoms in the cloud, but can, however, be explained within the Bose- 
Hubbard model 

The two site BH model predicts for the case of attractive interactions the exis- 
tence of a strongly correlated ground state for specific values of the parameters-i^. 
These ground states are far from being of mean-field type thus exhibiting inter- 
esting quantum properties, e.g. cat-state-like behavior. The existence of strongly 
correlated ground states of quantum systems has recently been linked to the exis- 
tence of instabilities of the semi-classical predictions in several different contexts: 
sonic analogues of black holes vortex nucleation in small atomic cloud s 
Bose-Einstein condensates (BEC) in rotating ring superlattices^^, or in the ground 
state of BEC in a double- well potential^-. 

The simplicity of the two-site BH model allows for an exact numerical solu- 
tion. However, it is always useful to have analytical insight which captures the 
essential physics sometimes hidden in the numerical diagonalization process of 
the Hamiltonian. In Ref. [ 191 . a mean-field state was proposed and the need to 
go beyond this approach was clearly established especially in the bifurcation re- 
gion where the cat-like states were identified. Later works along the same lines, 
but more focused on the dynamical properties, have contributed to the analysis of 
these systems ^^». Based on the limitations of the existing variational states, in this 
manuscript we propose an improved variational ansatz which is shown to yield 
an accurate description of the exact ground state of the system for a broad range 
of interaction strengths, including the strongly correlated regimes. This improved 
trial state is constructed by combining two states of mean-field type, thus also pro- 
viding an analytical representation of the ground state of the system, capturing its 
most representative features. 

The manuscript is organized as follows. In Sect. |2] we recall the definition of 
the two-site Bose-Hubbard model and introduce the tools to analyze the system. 
We also comment on some standard results obtained by exact diagonalization. Fol- 
lowing the steps of Ref. [ 19], in Sect.[3]we analyze the advantages and limitations 
of the mean-field approximation, and explore the possibilities of a variational state 
that describes the system beyond mean-field. In Sect. |4] we propose an improved 
variational state that can be used in the full range of the interaction strength and 
that incorporates the mean-field description, the incipient fragmentation before the 
bifurcation and the strongly correlated cat-like states whenever they are present. 

The main results and the conclusions are summarized in the last section. 



2 Theoretical description 



A good description of a system with particles that populate two weakly coupled 
states, which could represent the two sides (left and right) of a double-well, and 
with weak interaction between the particles that occupy the same state, is provided 
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by the Bose-Hubbard model, 

H = -e{a^i^ai- u^ur) ~ J{a\aR + aia^i^) + ^{a^ia\aLaL + a\a\aRaR) . (1) 

Where / describes the coupling between the two states, i.e. tunneling in the case 
of a double-well. Here, U characterizes the interaction between the particles and is 
taken to be the same in both sites. U > {U < 0) describes a repulsive (attractive) 
interaction. A small bias, < £ ^ /, is introduced to ensure the breaking of the 
left-right symmetry. Positive values of e promote the L state Q. 

A natural basis to study the system is the Fock basis, which is character- 
ized by the number of atoms in each of the two modes, \Nl,N]{). This basis, 
{|7V,0), |A'^— 1, 1), ...,\1,N - 1),\0,N)} spans an + 1 dimensional space, where 
N = Ni+Ni{ is the total number of particles. 

The action of the creation and annihilation operators on these states is de- 
fined in the following way: ci^i^Ni^Nr) = \/Nl + 1 \Nl + 1,Nr), and aL\NL,NR) — 
s/NI\Nl ~1,Nr). Therefore, 

In the two-mode approximation, a general N-body state can be written as 

\^) = Y,c,\k,N-k), (3) 

and the average number of atoms in each mode for a given state is A'^ = (f|a^a^ |f^), 

with P =L,R. The population imbalance of a state \ W) and its dispersion are de- 
fined as: 

z={W\Z\W); a, = ^{W\Z^\W)-{W\Z\W)^ . (4) 

with Z = {a[ai ~ aj^aR) /N. 

To characterize the degree of condensation of the system one can make use 
of the one-body density matrix, p. For a state |f ), we have p,, = (f |p,,|f ), 
with pjj = ajaj and j = L,R. The trace of p is normalized to the total number 
of atoms. The two normalized eigenvalues, i.e. eigenvalues divided by the total 
number of atoms A', of p are «i(2), with rii > «2- They fulfill ni +«2 = 1- The 
eigenvalue «, corresponds to the condensate fraction in the macro-occupied single- 
particle state I \ffi) which is the i-th eigenvector of the one-body density matrix. 
When the eigenvalues of the density matrix are strictly n\ = \ and «2 = 0, the 
system is fully condensed in a single-particle state (eigenvector of p). In this 
case, it is possible to express | f ) with a mean-field state constructed as — 
|V/i)®...®|v^i> = IV^i)^^. 



Note that from here on in our discussion, we will use the nomenclature of two sites or two 
wells when we refer to the two weakly coupled states that define our Bose-Hubbard model. 
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Fig. 1 (Color online) Energies of the lowest energy levels with respect to the ground state energy 
£0, as a function of the parameter A = NU /J, for A' = 50 (left panel) and A' = 500 (right panel). 
All energies are measured in units of J. 



2. 1 Exact spectral properties 

In this work we fix 7 = 1, which is equivalent to measuring the energy in units of /. 
We will vary the number of particles A' and the strength of the interaction, which 
will be considered always attractive, U <Q. The bias term, that can be related 
to possible small asymmetries of the external potential, will be taken very small: 
£//= 10-^. 

In principle, one can calculate the matrix elements of the Hamiltonian in the 
Fock-space basis and by diagonalization obtain the full spectrum of the system and 
the spectral decomposition of the eigenstate s 19.26,7,28,29 _ i\)s,n straightforward 
to calculate also the population imbalance and the degree of condensation of each 
state. 

In this paper we are interested in understanding the physical nature of the 
ground state of the system, which is for some parameter values quasi-degenerate 
with the first excited state. Therefore, we start by considering the lowest energy 
levels of the system as a function of A = NU /J, a parameter governing the behav- 
ior of the system. In Fig. [l] we report the energies of the first three excited states 
with respect to the ground state of the system as a function of A for two different 
numbers of particles, N = 50 and A'^ = 500, obtained by direct diagonalization—. 
For vanishing atom-atom interaction, A = 0, the energy gap for consecutive states 
is equal (except for the bias), and the gap is independent of the number of particles. 
As |A I increases the eigenvalues start to merge in pairs (the ground with the first 
excited, the second with the third, etc.) but due to both £ and /, they do not reach 
complete degeneracy. Moreover, the convergence of the merging process depends 
on the number of particles: for higher A' it occurs at smaller values of | A | reaching 
the value |A | = 2, when the number of particles tends to infinity. 

In Fig.|2](a) we plot the spectral decomposition of the ground and first excited 
states in the Fock space for different values of A, and for A'^ = 50. The plotted 
values \c]t\'^ give the probability that the state has k particles in the left well and 
N — k particles in the right one. Notice that if the spectral decomposition of the 
state is peaked at high values of k, it means that for this state most of the atoms 
are located on the left side of the double-well. For weak interactions, |A | < 2.6, 
the spectral decomposition of the ground and the first excited states are clearly 
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Fig. 2 (Color online) (a) Spectral decomposition (Iq^) in the Fock space of the ground (black 
solid line) and first excited (red dotted line) states for different values of A , with TV = 50. To help 
in the reading of the figure, instead of plotting the discretized values |qP we have generated a 
smooth curve by joining the different points, (b) Population imbalance z (black solid line) and 
its dispersion 0%, Eqs. (black dashed line) of the exact ground state as a function of A. The 
semi-classical predictions of the imbalance (red dotted line) and its dispersion (red dot-dashed) 
are also plotted (c) Solid and dotted lines depict the condensed fractions n\ and n2 of the 
one-body density matrix of the exact ground state as a function of A. In all cases A' = 50. 



different, (as were also the energies in Fig.[T])- For stronger interactions, —3.2 < 
A < —2.6, the two states become very close in energy (Fig. [U left panel) and 
their spectral decompositions are very similar. However, one should notice 
that the ground state is symmetric, c^/2-i-<: = <^A'/2-*:' ^nd the first excited one is 
antisymmetric, c^/2+*: = ^<^w/2-j(:- In this region the ground state is a strongly 
correlated cat-like state0 as its spectral decomposition has two clear peaks. 

Finally, for |A | > 3.2, the two states become again clearly different: the ground 
state is peaked at a high value of k, with a large amount of atoms in the left well, 
while the first excited has its peak at a low value of k. Note that the energies of 
these states are very close to each other. 

A useful characterization of the ground state is provided by the population 
imbalance z- As shown in Fig.|2](b), it remains zero up to a certain value of |A| 

3.25 for N = 50), approaches 1 as |A | increases further. The figure also shows 
a^, which starts from small values associated to a relatively narrow binomial dis- 
tribution. It increases in the range where the strongly cat-like state is present, and 
finally decreases abruptly when | A | increases further and the ground state popu- 
lates massively the L state. Thus z — !> 1 and — !■ for |A | — 00. 

The degree of condensation of the ground state, |(/)gs), is determined by the 
eigenvalues n\ and «2 of the one-body density matrix, which are plotted in Fig.|2] 
(c). These condensate fractions measure the macroscopic occupations of the single- 
particle states and Iv'i), eigenfunctions of (0gs|p|(/)gs). The regions where 
these values are not close to 1 and 0, signal the occurrence of fragmentation 
of the ground state and the impossibility to describe the system by means of a 
mean-field state. In the region, —2 < A < 0, «i is rather close to 1 {n\ ^ 0.99), 

^ Strictly speaking, the purest cat-state would correspond to the state l/\/2(|A',0) + |0,A')). 
The states we refer to as cat-like states are sometimes called kitten states with a certain degree 
of 'catness' 
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and the macroscopically occupied state is given by | y/i) = + \R))/^/2, with 
\L{R)) = a^^jlO). However, as we will discuss later in Fig. [3] this slight frag- 
mentation produces noticeable differences in the spectral decomposition of the 
mean-field state build with the state | ) and the exact ground state. The fragmen- 
tation is particularly important for —2.5 > A > —3.5, which is roughly the same 
interval where the cat-like structure takes place. However, the macro-occupied 
state I y/i) remains equal to the one previously discussed. The correlations beyond 
mean-field affect the degree of condensation, but not the single state that is mainly 
occupied. This is because the ground state remains almost symmetric (except for 
the bias) in the Fock space (see Fig. |2](a) andO, i.e. with z almost zero. This is 
reflected in the symmetric character of the one -body density matrix, which in turn 
implies that is the normalized symmetric combination of \L) and \R). For 
further increasing |A|, the system becomes again condensed: ni — > 1. The slight 
energy difference introduced by the bias term, which energetically promotes the 
\L) over the \R) state, drives the system to li/i) — !■ \L). 

The precise value of the bias term has been shown to determine, for a fixed A'^, 
the size of the cat-like -region. Exploring the interplay between the bias term and 
the hopping strength, /, a good estimate of the precise value of A where the bias 
dominates is given in Ref. [ 30]. For larger number of particles, A^, the bias term 
becomes dominant at lower values of |A | as its effect is proportional to A'^, and 
therefore the cat-like-region becomes narrower at values of A closer to the critical 
classical value A = 2. 



3 Variational state for the ground state 

3.1 Mean-field ansatz 

A reliable mean-field statei^ can be constructed using a general single-particle 
state |(/))sp = a|L) + j3|/?), with |a|^ + |j3|^ = l, and considering all the particles 
to be in this single-particle state: 



The expectation value of the Hamiltonian for this state is, 

E[a,a\p,P*) = {^\H\^)N = -eN{aa* -l5l5*)-JN{a*l5 + al5*) 

+|Ar(7V-i)(|a|4 + |/5|4). (6) 

The minimization of the energy with respect to the variational parameters, together 
with the normalization condition |ap + |j3p = l, yields the following equation 

2£N-Jn{^^^^^ +UN{N-\)m^-\a\^)=0. (7) 

The possible solutions of the previous equation will be of the type (a,±j3) with 
both a and j3 positive real numbers. Explicit simple analytic solutions to the previ- 
ous equation can be obtained by neglecting the bias term. Therefore, taking e = 
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and introducing A =A{N ~ l)/N, one gets the following set of solutions^: 



2^1/22 



(8) 



that give rise to the multi-particle states: 



\^t)N = -i^_[cCial±l5ialY\0), 



(9) 



with « = 0, + , — . Note that the solutions a± and j3± only exist when |A | > 2. The 
expectation value of the energy in these states, with £ = 0, is : 



= {^^\H\^^) = ^N{N-1)tJN, 



(10) 



^ (0^ =E^^ ieme) = \n{N-\) ^(1 t2) . (11) 



As C/ < and / = 1, the states \ (^^)n have a lower mean energy than the |(/>~)a? in 
all cases. To find the lowest energy, we study the difference between E^ and £^ 
(notice that = E^), 



E+ -E^ 





"1 - 




NJ 













(12) 



Thus for |A| < 2, the lowest energy state is \^q)n, and for A < —2, both states 
\^^)n and |0j^)Ar have the same minimum mean energy. 



3.2 Variational ansatz beyond mean-field 

In the case that E^ are the smallest mean-field energies, i.e., for \A\ > 2, one can 
propose an alternative ansatz for the multi-particle many-body state that goes 
beyond the mean-field approach and tries to incorporate the cat-like structure!; 



'cat//V 




The expectation value of the Hamiltonian for this many-body state. 



(0cat|/^|^>cat)iV = 



NJ 
4A 



4-A' 



l + (2/|A|r 



3A' 



(13) 



(14) 



is smaller than E^ . This state is a linear combination of two non-orthogonal mean- 
field states having the same energy expectation value (if the bias is not taken into 
account), but two different spectral decompositions. It is precisely the fact that 
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Fig. 3 (Color online) Spectral decomposition (q) in Fock space of the ground state of the sys- 
tem, for different values of A , computed by exact diagonalization (black solid line) and com- 
pared to the spectral decomposition of the mean-field functions |(^q")w for |A| < 2 (red dashed 
line) and |<^+)w for |A| > 2 (blue dashed line). For |A| > 2, we also show the results for the 
variational cat-state |(|>cat) (green dotted Une). In all cases A' = 50. 



they are not orthogonal that allows the mean energy value in the state | (/>cat) to be 
smaller than E'^ = E^. 

In Fig. [3] we show the Fock space decomposition (c^) for different values of A 
of the ground state of the system computed by exact diagonalization of the many- 
body Hamiltonian, Eq. ([T|. We compare these coefficients with the ones provided 
by the mean-field state \^q)n for \A\ < 2, and \^^)n for |A| > 2. In this last case 
we also plot the results for the variational state \(j)cat)N given in Eq. ( 1131) . In all 
cases A'^ = 50. Note that for this number of particles A and A are very similar 
and therefore the critical value of A where the mean-field states |^±)/v appear is 
A - A = -2. 

For |A I < 2, the best mean-field representation of the ground state corresponds 
to \^q)n- The coefficients Ck follow a binomial distribution, symmetric around 
k — N/2. This mean-field state gives a good qualitative description of the system 
in this range, however it coincides with the exact solution only for A = 0. In Fig.|3] 
one can appreciate that the distribution of the exact ground state is slightly broader 
and the differences increase with |A | (recall that the fragmentation in this region 
was very small). The energy difference between {^^H\^q) and the exact ground 



state energy, E^^, relative to E^^ is shown in Fig. |4] (a). This relative difference 
increases with A and is zero only for A = 0. Another measure of the capability of 
the mean-field state to describe the exact ground state is provided by the overlap of 
the trial state with the exact ground state. This overlap, ((/)gs | (/(o^) is plotted in Fig.|4] 



^ Note the only difference with the state used in Ref. l ITgll is due to our state being normalized 
to 1. 
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Fig. 4 (Color online) (a) Relative difference with respect to the exact ground state energy of the 
expectation value of the Hamiltonian with different variational many-body states: \(j>g} (violet 
dotted line), (blue triangles), \tj)cat) (Eq.[T3} (green dot-dashed line), |'f^ai-)min (red dashed 
line) and |'Ivar)max, (black solid hne). (b) Overlap of the different states discussed in the text 
with the exact ground state of the system as a function of A . | (^q") (violet dotted line), | (j)^) (blue 
triangles), (black crosses), \^cRt) (Eg.fTst (green dot-dashed line), |'f^ar)min (red dashed line) 
and |'f^ar)max (black solid line). 



(b). The overlap is 1 only for A = 0. The differences in the spectral decomposition 
in the Fock space reflect in an overlap smaller than 1 when A increases. Clearly 
at A = —2 the overlap between ((/)gs|(/)Q^), decreases quickly and tends to zero for 
large values of |A|. 

In the region — 3.2 < A < —2, the minimum energy mean-field solutions, | (j)^) 
and \(pl), provide the same energy expectation value when the bias is not taken 
into account. The Fock decomposition of is plotted in Fig. [3] The distribution 
for would be symmetric with respect to k — N/2. The differences of the ex- 
pectation energies, of these two states ) and |<^^ )), with respect to the ground 
state energy are rather small, not only in the cat-like state region but also for larger 
values of |A |, where the difference tends to zero. On the contrary, the behavior of 
the overlap of these two states with the ground state is rather different. In the cat- 
state region, both overlaps are rather similar. The reason is that \ (j>^) overlaps with 
the right part of the cat-state (in the Fock space) and |^^) overlaps with the left 
part of the cat-state. As |A | is increased and the cat-state disappears, the presence 
of the bias term in the Hamiltonian ensures the breaking of the left/right symme- 
try by energetically promoting the \L) state. Thus, the system selects |^^) as the 
ground state, and therefore its overlap with the exact ground state tends to 1, while 
the overlap of ((/)gs|(/)^) tends to zero. 

The cat-state structure can be reproduced by defining as trial state the linear 
combination of \(j>^) and Eq. ( 113b , as one can see by looking at the spectral 
decomposition of this state shown also in Fig. [3j Obviously, when | A | increases, 
and the ground state is preferentially located in one of the wells, this state | (j)cat) 
does not give anymore a good reproduction of the Fock space decomposition of 
the ground state. If one looks at the overlap (^gsl'/'cat), this variational state clearly 
improves the overlap with the ground state in the cat-like region, but when |A | 
increases, the overlap tends to a constant 1 / \/2. The behavior of the energy can 
be observed in Fig.|4](a). We can see that there is an improvement in the cat-state 
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region when using |^cat)- However, when \A \ increases all three functions |^cat), 
|(/)^) and \ become degenerate in energy with the exact ground state. 



4 Improved global variational ansatz 

We propose a variational ansatz that is valid independent of the strength of the 
interaction including at the same time the possibility of a mean-field and the exis- 
tence of a cat-state. This state is a combination of two different mean-field states: 



|fvar) =A|(/>)i 

B 



/Nl 



fB|(/))2 = 

/3a| + aa[ 



N 

|0) 



aa]+^ai |0) 



(15) 



The variational parameters a, j5, A and B are taken real. Note that the two mean- 
field states are not necessarily orthogonal and therefore the normalization condi- 
tions are imposed in the following way: 



a 



15^ = 1, A^+B^+2{2al5fAB=l. 



(16) 



Let us discuss the differences of the ansatz in Eq. dlSI ) with the states studied in the 
previous section. Here if A or S are zero, the state reduces to a mean-field state of 
the type considered before. On the other hand, if one constructs the combination 
of the two mean-field states and allows for a new minimization of the variational 
parameters, a noticeable improvement of the state is obtained. The expectation 
value of the Hamiltonian with this ansatz is given by 



E 
JN 



-lap 



A 



l-4a2j32) 



1 



{l5^-a^){A^-B^)- 



2aj3 
A 1 
'4 



{A^+B^ 



lap ■ 



(17) 



To determine the parameters of the variational state we follow two different 
criteria. The first consists in performing a numerical minimization of the expecta- 
tion value of the energy, Eq. ( 117b . The many-body state thus computed is named 
|fvar)min- In the second procedure, which can be pursued only when we already 
have a numerical solution of the exact ground state, we determine the coefficients 
by maximizing the overlap of the variational state with the exact ground state, 
giving the state |fvar>max- 

The first procedure, which does not require the previous numerical solution of 
the ground state, produces by construction the closest energy to the exact ground 
state energy within the form of Eq. dlSb . As will be discussed in the following, the 
second criteria although requiring the previous numerical solution of the ground 
state, produces in all cases an extremely close agreement with the ground state 
from the energetic point of view, while also improving the overlap with the numer- 
ically computed ground state. Thus, for certain applications where an analytical 
rendition of the state is preferable, our variational proposajf] should be very useful. 



In the sense that (f^a 



ir)max provides also an upperbound to the ground state energy 
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Fig. 5 (Color online) Values of a, j8 (a) and of A, B (b) obtained in the improved global 
variational approach of Eq. il5\ by overlap maximization (black solid line) and by energy min- 
imization (red dashed line), and the ones obtained in the variational approach of Eq. <13t . for 
which A = B (green dot-dashed line). 
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Fig. 6 (Color online) Spectral decomposition of the exact ground state (black solid line), and 
the states obtained using the ansatz defined in Eq. i fTsl l when its overlap with the exact ground 
state is maximized (green dashed line) or when its energy is minimized (red dotted line). 



The values of a, j3, A, and B obtained in both cases are reported in Fig.|5] For 
A = 0, we have a = P = 1/ \/2 and A = S = 1 /2, recovering the function | (p^) 
that was the exact solution. This is the only case where our improved variational 
state coincides with the one proposed in Ref. [ fl9ll . Obviously in this case both 
conditions: minimum energy and maximum overlap, provide the same solution. 
Note that in this case, the overlap between the two components of the generalized 
variational ansatz and |(/))2) is maximum, i.e. the two components coincide. 

When I A I is increased, — 2 < A < 0, a and /3 become different while A and B 
remain equal but different from 1/2. The improved variational state incorporates 
correlations beyond mean-field, and the overlaps ((/)gs|'fvar)min and ((/>gs|'f^ar)max, 
[see Fig. |4](b)] clearly improve with respect to Also the difference in the 

expectation values of the energies (f^arl^^l'fvar)min and (f^ai |^^|^ar)max relative 
to the ground state energy become smaller as shown in Fig.|4](a). 
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Fig. 7 (Color online) (a) Population imbalance z and (b) its dispersion a- as a function of A 
for the exact calculation (black solid line), the improved global variational approach in Eq. (Tsjl 
by overlap maximization (red dashed line) and by energy minimization (violet dotted line) The 
blue dashed line is the semi-classical prediction and the green dot-dashed line corresponds to 
the |(^>cat)A' of Eq. i fTst . The number of particles isN = 50. 



In this region (|A | < 2) the differences between the observables corresponding 
to these two variational states associated with the maximum overlap or with the 
minimum energy criteria are rather small. The state that minimizes the energy pro- 
vides slightly better energies, however this difference is not significant in Fig. |4] 
(a). Correspondingly the state that maximizes the overlap provides overlaps with 
the ground state closer to unity. However these differences are also not appreciable 
in Fig.|4](b). The Pock decomposition of these two variational states |f^ar)min and 
|'fvar)max Compared with the one of the ground state are shown in Fig. |6] for dif- 
ferent values of A and = 50. One can observe a clear improvement of the Pock 
decomposition with respect to the mean-field state in this range of the interaction 
|A|<2. 

Interestingly, the proposed state captures well the correlations beyond mean- 
field existing in the ground state of the problem before the classical bifurcation. 
These correlations, as discussed above and shown in Fig.|2](a), produce very small 
effects on the condensate fractions but become clearly visible when looking at 
the spectral decomposition of the ground state, see Fig. [3] or the dispersion of 
the population imbalance which is no longer corresponding to a simple binomial 
distribution, see Fig.[2](b). 

Once we cross the classical bifurcation, A < —2, the spectral decomposition 
of the ground state broadens and at A —2.2 becomes quickly two-peaked. This 
region where the ground state has two maxima is what we refer to as the cat- 
state region. The main objective of the variational ansatz introduced in [ 19], and 
discussed in Sect. 13.11 is to describe the ground state properties in this region. 
The results with the improved global variational ansatz of Eq. ( 115b are shown in 
Figs. |4] |6] and [7] Unlike in the region before the bifurcation, here the two crite- 
ria used to compute the variational parameters provide fairly different results in 
some cases. The computed energy of the state is very close with both criteria, but 
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Fig. 8 (Color online) Largest eigenvalue of the one body density matrix, m, as a function of 
A for the different many-body states discussed in the text: exact calculation (black solid line), 
|^var)max (red dashed line), |^var)min (violet dotted line), and the cat-state |(^(.at) of Ref. [ [I^ 
(green dot-dashed line). The inset shows the region before the bifurcation. 



its overlap with the ground state is different depending on the criteria used, see 
Fig. m This is a consequence of the clear differences seen in the spectral decom- 
position. Fig. |6] The variational solution obtained by minimizing the energy is 
seen to depart from the exact solution in the region —3.5 <A < —3, predicting an 
earlier transition to the 'self-trapped' domain, see Figs. |5] and |7] 

The criterion of maximizing the overlap implies the ansatz to follow much 
closer some of the explored ground state properties of the system as shown in 
Fig. |7] where the agreement with the exact calculation both for the population 
imbalance and its dispersion is extremely good in the considered domain. Thus, 
for these parameter values, obtaining a faithful representation of the ground state 
with the form proposed in Eq. ( 115b would require the prior numerical solution of 
the exact ground state. 



4. 1 Fragmentation of the ground state 

To complete the characterization of the proposed states we also study the frag- 
mentation of the ground state of the system. To this end, we calculate the one- 
body density matrix and look at its larger eigenvalue. If the largest eigenvalue (rii) 
is significantly smaller than unity, we have fragmentation and the system is not 
condensed in one single state. It also indicates the impossibility to describe the 
system by a mean-field state and therefore reveals the existence of correlations 
beyond mean-field. The largest eigenvalue of the one-body density matrix associ- 
ated to the different states discussed in this work is reported in Fig.|8]for different 
values of I A I . The two mean-field states, (| (/)q^) , | (/)^)), are not plotted as they have 
this eigenvalue equal to unity independently of A . The exact ground state gives 
rise to an rii very close to unity, in the region |A | < 2. However, the eigenvalue is 
strictly one only for A = 0, and is actually a smooth decreasing function of | A | . 
It decreases faster in the cat-like region reaching a minimum (ni ^ 0.8) (maximal 
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fragmentation) around A = —3.2. For larger values of \A \ it grows again reaching 
the value 1 as the system condenses in the left well due to the bias. 

The Hi associated to the variational state, |(/)cat), which exists for \A \ > 2, starts 
from wi = 1 and decreases with increasing |A | reproducing rather well the exact 
111. Contrary to what happens with the exact tii, it continues decreasing and in- 
creasing the fragmentation failing to reproduce the region dominated by the bias. 
Finally, the variational many-body states proposed in the present paper, |'S'var)min 
and |'S'var)max reproduce very well the exact «i in the region before the bifurca- 
tion, where the system is slightly fragmented, see the inset in Fig. [H This small 
fragmentation, as discussed above, indicates the presence of some correlations 
beyond the mean-field already in this region. In the cat-state region, the |<S'var)max 
also reproduces the exact ni, see Fig. [8] 

5 Conclusions 

The variational analytical approach to the two-site Bose-Hubbard model gives a 
useful insight into the physical nature of the ground state of this apparently simple 
system that however shows a very rich phenomenology when the interaction or the 
number of particles change. We have carefully studied the limitations of the mean- 
field description strongly linked to the presence of fragmentation of the condensate 
and quantum fluctuations. The proposed variational state is able to describe rather 
well the exact state and reproduces the energy, the one-body density matrix and 
thus the fragmentation of the state which are the main magnitudes that we have 
used to characterize the ground state. 

We have also compared the spectral decomposition of the exact ground state 
with the proposed state obtaining good agreement. The many-body states | <S'var)mini 
whose parameters are obtained by minimizing the energy, can be used for any 
number of particles. This state, incorporates for all As, quantum correlations be- 
yond the mean-field and reproduces very well the fragmentation induced by these 
correlations which become larger in the cat-state region. 
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